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jrt Abstract. In this paper a length-conserving numerical scheme for a nonlinear fourth order sys- 

^^H tern of partial differential algebraic equations arising in technical textile industry is studied. Ap- 

| .^ l plying a semidiscretization in time, the resulting sequence of nonlinear elliptic systems with alge- 

braic constraint is reformulated as constrained optimization problems in a Hilbert space setting 
that admit a solution at each time level. Stability and convergence of the scheme are proved. 
The numerical realization is performed by projected gradient methods on finite element spaces 
which determine the computational effort and approximation quality of the algorithm. Simula- 
tion results are presented and discussed in view of the application of an elastic inextensible fiber 
motion. 
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1. Introduction 



^ 



^ The numerical simulation and optimization of the dynamics of thin long elastic fibers are of 

^^ great importance in the technical textile industry, for example in production processes of yarns or 

rij non- woven materials [HI [TT] . In the slender-body theory a fiber can be asymptotically described 

fy-N by a parameterized, time-dependent curve that might represent the center-line. Then the fiber 

"^ dynamics can be modeled by a nonlinear partial differential algebraic system (PDAEs) [T3]- The 

(—» K partial differential equation of fourth order for the curve has a wavelike character with elliptic 

rvq regularization. It can be considered as a reformulation of the Kirchhoff-Love equations for an 

T— I elastic beam [121 E]- The nonlinearity enters the system by the attachment of an algebraic constraint 

IJ that incorporates inextensibility into the model and turns the inner traction to an unknown of the 

• ^H system, i.e. Lagrange multiplier. Depending on the application, also nonlinear or even stochastic 

/\f source terms might play a role, see [iS] and Figure |1.1| The studies of fiber lay-down processes 

and longtime behavior require a fast, accurate numerical treatment, cf. [TT1[51IH]. So far, the used 
approaches were mainly addressed to high-speed performance without any theoretical results on 
convergence or length conservation. 

This work aims at the development of an appropriate numerical scheme for the PDAEs with 
focus on analytical and computational aspects. We propose a semi-discretization in time. Following 
the concept of |10| and employing a horizontal line method, we replace the transient problem by 
a sequence of elliptic systems that are handled in their weak formulation in terms of the Lagrange 
formalism. The algebraic constraint is incorporated in the definition of the optimization domain 
such that we study the solvability of a constrained minimization problem in a Hilbert space setting 
[SJ [TB]. We prove the existence of the minimizer and of the Lagrange multiplier on each time 
level. Stability estimates on the discrete solution and the Lagrange multiplier result then in the 
convergence of the numerical scheme. In addition, we derive an error bound on the fiber elongation. 
Numerically, we solve the optimization problems in finite element spaces by applying a projected 
gradient method with Armijo step size rule. Thereby, the finite dimensional approximation of the 
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Figure 1.1. Application: melt-spinning process of non- woven materials. From 
left to right: Turbulent air flow in process (photo by industrial partner), mean 
velocity flow field, turbulence effects on immersed fiber modeled by stochastic forces 



(source terms) in fiber PDAEs (2.1 1, 14 



constraint and thus the choice of the projection mapping determine the efficiency and accuracy of 
the scheme. 

The paper is structured as follows. We start with the presentation of the model equations in 
Section [2] After the introduction of the numerical scheme we deal with its theoretical analysis - 
regarding existence, stability and convergence - in Section [Sj The numerical realization is discussed 
in view of computational effort and approximation quality in Section |4J The simulation results 
are in accordance with the analytical investigations and satisfy the predicted error estimates. We 
conclude with a summary and an outlook. 



2. Model 
Due to its slender geometry, a fiber can be asymptotically represented by its arclength-parame- 



terized time-dependent center-line r : il — ili x fix 



according to the special Cosserat theory 



0, where fla '■= (0,a), a g (0,oo) with fiber length I and end time T. Since extension and shear 
are here negligibly small in comparison to bending, the dynamics of an homogenous fiber can be 
described by a nonlinear partial differential equation of fourth order with an algebraic constraint 
for the conservation of length [T31 



ujdttr{s,t) = ds{X{s,t)dsr{s,t)) - bds 
|9.r(s,t)p = l 



r(s,i)+f(r(-),s,0 



(2.1a) 
(2.1b) 



where a; > denotes the line weight. The dynamics is caused by acting outer and inner forces (New- 
ton's law). The outer force densities f might come for example from gravity and/or aerodynamics 
regarding the application. The inner force densities stem from traction A and bending with bending 
stiffness 6 > 0. The inner traction A : fl — >■ K might be interpreted as Lagrange multiplier to the 
algebraic constraint (2.1b) that is expressed in the Euclidian norm | • |. In addition, torsion could 
be included, yielding an extra term K{dsY x dsss^) in (2.1a) with dsK — 0, [TT]. Its neglect {k = 0) is 



here justified by the set-up that is characterized by a free fiber ending. In particular, we consider a 
fiber fixed at one ending [s — I) that is freely swinging. This set-up is clearly a simplification to real 
applications in technical textile industry, but it still contains the major mathematical difficulty, i.e. 
the partial differential-algebraic structure of the model equations. Then, Dirichlet and Neumann 
boundary conditions for the fixed and stress-free fiber ending respectively as well as appropriate 



initial conditions close (2.1) to an initial boundary value problem 



r(/,i)=0, 
a,,r(0,i) = O, 

r(s,0) = (;-s)eg. 



atr(s,0) = O. 



A(0,t) = 0, 



(2.2) 



with unit vector Gg prescribing the direction of gravity. 
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The fiber model (2.1a) is a wavelike system with elliptic regularization. It is a reformulation 



of the KirchhofF-Love equations for an elastic beam [T^]. The mathematical challenge lies in the 



treatment of the algebraic constraint of inextensibility (2.1b). This paper aims at the derivation and 
investigation of a numerical scheme for the partial differential algebraic system. For simplicity, we 
restrict here to sufficiently smooth outer forces that are independent of the fiber position, f : fi — >■ K'^, 
like gravity. However, note that in melt-spinning processes of non-woven materials, even stochastic 
forces from turbulent air flows become important [TSIIII], cf. Figure [lT| For a study on stochastic 
extensible beam equations see e.g. [HllSj- 

3. Theoretical analysis 
In this section we derive an implicit semi-discretization for the partial differential algebraic fiber 



system (2.1 1. Following the concept of [10] and employing a horizontal line method (Rothe method) 
in time, we replace the transient problem by a sequence of elliptic systems that are handled in 
their weak formulation in terms of the Lagrange formalism. Hereby, we incorporate the algebraic 
constraint in the definition of the optimization domain such that we study the solvability of a 
constrained minimization problem. In particular, we show the existence of the minimizer and of 
the Lagrange multiplier on each time level. Stability estimates on the discrete solution and the 
Lagrange multiplier result then in a convergence proof for the numerical scheme. 

3.1. Semi-discretization. Let T e (0,oo) be given. We divide the time interval [0,r] into N 
subintervals by introducing the temporal mesh {t^ | fc = 0, ..., N} where tk = kr is prescribed by the 
time step r = T/N. Instead of the taken equidistant discretization one can also think of a varying 
time step Tk = tk — ife-i and define the maximal subinterval length t = maxfc=i_...^7v t^. Then the 
partition is assumed to satisfy r — >■ as TV — >■ oo. 

Using an implicit Euler scheme, we discretize the fiber system as 

'^^ "^ ds{Xk+idsrk+i)-bdssssrk+i+ik+i (3.1a) 



IdsTk+ll^^l, (3.1b) 

with tk+i — f (ife+i), k = 1, . . . ^ N — 1, Tq = (I — s)eg and ri = Tq. The implicit time discretization 
requires consequently the recursive solving of nonlinear elliptic systems in one space dimension. The 



approximate solution to (2.1) is then given by the linear interpolation (r'^, A"^), i.e. 



r^{s,t) = - — !-(rfe -rfc-i) +rfe_i, seili, t e {tk-i,tk], k = l,...,N, 

T 

and correspondingly for A'^, where Aq = Ai = 0. For functions defined on [0, T], in turn, a subindex 
A: e {0, . . . , N} corresponds to the value of the function at time tk- 



The discretized system (3.1 1 can be identified as Euler-Lagrange equations corresponding to the 
appropriate Lagrange functional such that we study its solvability as variational problem in the 
following. 

3.2. Solvability of discretized system. The norm of a Banach space B we denote by |1 • ||b and 
the dual pairing with its dual space B' by b{'i')b>- If ^ even is a Hilbert space, then its inner 
product we denote by (•,-)s- By L^(S';R"), S C 'R'^ Lebesgue measurable, n,d & N, we denote 
the Hilbert space of (equivalence classes of) square integrable functions on S w.r.t. the Lebesgue 
measure taking values in M". The space C°(S'; B) of continuous functions on compact S with values 
in B we consider to be equipped with the norm of uniform convergence. We use the notation 
W"''P{U;W), U CR"^ open, m e [0,oo), p e [l,cx)], for Sobolev spaces as in [I]. In case of p = 2, 
the Hilbert space I4^™'^([/;K") is abbreviated by H"^{U;R"'). In the case n = 1 we suppress the 
range of function spaces. In particular, we introduce the notation 

H^^{na;W") := {v e if™(ri,;M") | 9Xa) = for all a e No, a + 1/2 < m}, Qa = (0,a), 

a e (0,oo). Of course, HQ\{ila;M.^'') equipped with the norm of i/™(51a;M") is a Hilbert space. 
Its dual space (iJ^„(r2a;M"))' we denote by iJ-™(fia;K"). Recall that iJ™(f^a;M") is embedded 
continuously and compactly in the Holder spaces C''''*'([0, a];IR") for m > 1/2 + k + j, k E Nq, 
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< 7 < 1, see e.g. [T]. We always, via the Riesz representation theorem, identify spaces of square 
integrable functions with their dual space and consider an embedding of Sobolev spaces in the sense 
of Gelfand triples with the space of square integrable functions as central space. 
Then, we define the affine linear fiber space 



V:={vevD + HliiQi-.R^) I ^D G H\ni;R^), vd{1) = 0, 9,vd(0 = -eg}. (3.2) 

^ = and express it in terms of 

(3.3) 



We consider the constraint of inextensibility (3.1b) as dtldsr^^i 

efe+i : V -^ Hliini), efc+i(v) = 2ds{v - Vk) ■ 9,rfe = 



Moreover, we deduce the cost functional Jk+i '■ V 



^fe+i(v) = w||tD^,_^iv| 



L2(n, 



+ 6||a,^v||^2^j^^^ -2(ffe+i,v)i2(o,) 



(3.4) 



with the second temporal difference D^ , -^v = (v — 2rfc + rfc_i)/r^ by applying variational calculus 
on p.lal) for /c = l,...,iV- 1. 



Lagrange formalism. For fc = 1, ..., A^ — 1, let Jk+i be the cost functional of (3.4) and Ck+i the 
constraint functional of (3.3). Define the Lagrange functional L^+i '■ V x H^^(ili) — > M 6?/ 

Lfc+i(v,A) = Jfc+i(v) + ffi^(o,)(efc+i(v),A)^-i(j^^). 



Then, the minimizer of the Lagrange functional is a weak solution of the fiber system (3.1 ). 



The minimizer of the Lagrange functional satisfies the adjoint problem (3.5) for all test functions 
ri e H-\ni) and cj) e Hli{ni,R^), i.e. 

dxLk+i{^r,X)[r]] = = Hi,,(o,)(efe+i(v),?7)H-i(o,) (3-5a) 

VvLfc+i(v,A)[</.]=O = 4+i(v)[0] + ^i^^(o,)(e',+i(v)[(/.],A)^_,(^^, (3.5b) 

= 2(w (Dfe+iv, ct))L2{n,) + b {dssv, dss(p)L^(n,) - (ffc+i, 0)l2(o,) 

Presupposing sufficient regularity of the Lagrange multiplier A, the duality pairing ^i (o,)(-, ■) n-'^in ) 
coincides with (•, ■)L^(^n,) in the sense of a Gelfand triple. This yields the Euler-Lagrange equations 
corresponding to the fiber system. Hence, the weak solvability of the elliptic fiber system (3.1 ) can 
be formulated as 



Constrained minimization problem 

Minimize Jk+i over the domain Kk+i :— {"v £V \ Ck+i^v) = and \dsVk ~ 9sv| < r }. (3.6) 

Lemma 1 (Properties of cost functional). For k = 1, ..., N — 1, the cost functional Jk+i : V —>■ M. 
defined in (3.4) is strictly convex, coercive and weakly lower semi- continuous. The minimization 
domain K^+i is closed and convex and, in particular, weakly closed. 

Proof. Here and throughout the following proofs where is no danger of confusion, we suppress the 
indices indicating the time levels for a simpler notation. Let u, v e y, u 7^ v, /x e (0, 1). Then, the 
strict convexity of J is concluded from 



Ai J(u) + (1 - m) J(v) - J(/iu + (1 - /i)v) 
= (M-M')(^!kD2u-rD2v||i.(,,,) 



b\\dssn-dssA\^(ni)) > 



since w, 6 > 0. 

Due to the assumed boundary conditions a Poincare inequality holds and we obtain 



w||rD2v||^2(fi,) 



b\\dssm2^n,)>A\\v\\l,i^n^^-B 
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for some < A,B < oo. Hence 

J(v) - ^||rD2v||i.(f,^) + bWdssAh^n,) - 2{i,^)mn,) > ^l|v||l,.(j,,) - 2|(f, v)i2(o,)| - B 
> h\\H^(ni) (^l|v||//2(o,) - 2||f ||i2(o,)) - B. 
Thus, J(v) -^ oo, if ||v||//2(Qj) —T' oo for fixed f e L'^{ili), i.e., J is coercive. 

Let (v„)„gN be a sequence in V that converges weakly to v G 1^ in H , i.e., v„ ^ v for 

n —>■ oo. Then, in particular, v„ ^ v and 9ssV„ ^ 9ssV for n — > oo. Since the norm is lower 
semi-continuous w.r.t. weak convergence and the inner product with f e L-{Cli) is continuous w.r.t. 
weak convergence, we obtain J(v) < lim„_>oo inf J(v„), i.e., J is weakly lower semi-continuous. 

The convexity of K results from the afhne linearity of e and the convexity of the set {v e 
H^{Qi]R^) I la^rfc -9^v| < r^}. Let u,veK,fie [0,1], then ^u+(l-^)v e X, since 

e(/iu + (1 — /i)v) = A*e(u) + (1 — A*)s(v) = and 
l^^rfe - dsi^iu + (1 - /i)v)| < ^IdsTk - (9su| + (1 - M)|5si"fc - 9sv| < r^ 

Since V is closed, e is continuous and the set {v G H^{ni;R^) \ \dsrk — dgv] < r^} is closed, also 
K is closed. This, together with convexity, implies that K is also weakly closed. D 



Theorem 2 (Existence and uniqueness of minimizer). The constrained minimization problem (3.61 
has a unique solution Vk+i G K^+i on every time level. 

Proof. According to Lemma [lithe cost functional J on every time level is coercive and weakly lower 
semi-continuous and the domain K is convex and weakly closed. These are the requirements for 
a general existence and uniqueness result for constrained minimization problems, see e.g. [TB]. We 
state the proof here for completeness. 

Choose a minimizing sequence (v„)„gNj v„ G K, with J(v„) — > infveK ^(v) for n -^ oo. Then 
— oo < infvGK t/(v) < oo and (v„)„gN is bounded in view of the coercivity of J. Hence, there exists 

a subset N (ZH and r ^ V such that v„ ^^ r for N 3 n — > oo. Since K is weakly closed, r ^ K. 
The weak lower semi-continuity of J implies 

J(r) < inf J(v„), 

whence r is a minimizer. 

Since K is convex, the strict convexity of J on y implies the uniqueness of the minimizer. 
Assume u, v G /-C to be two minimizers, u 7^ v with J(u) = J(v). Then J(/iu + (1 + /i)v) < 
/^J(u) + (1 — fJ')J{^) = J{u) for fi G (0, 1). Since ^u + (1 + /.j)v G K for // G (0, 1), this contradicts 
the assumption. D 

Note that the uniqueness of the minimizer is meaningless for the solvability statement of the fiber 
system, since the unique minimizer need not necessarily be the only solution in view of possibly 
existing saddle points. 

Since 

\dsrk+i\ < IdsTk+i ~ d,rk\ + la^r^] < t^ + \d,rk - 5,rfc_i| + |9,rfc_i| < . . . < Nt^ + \e^\ = 1 + Tr 

and efc+i(rfe+i) = for all fc = 1, ..., A^ — 1, the following lemma follows immediately by the Cauchy- 
Schwarz inequality. 

Lemma 3. The relation 1 < |9srfc(s)| < |9srfc+i(s)| < 1 + Tt holds for all s G fli, k — Q, ..., A^ — 1. 

Proposition 4 (Surjectivity of linearized constraint functional). For k = 1, ...,N — 1 the linearized 
constraint functional ej.,]^ G L(ifg ;(r2/; M"^); iJg ;(r2j)) is surjective. 

Proof. We have 
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Let ip £ Hqi{Qi) be arbitrary. Set 



^i^) = - I wr^r^du, sen,. (3.7) 



2 1 duffel 

Note that dsV^^ip € C"([0, Z]). Hence, together with Lemma pi we obtain <p e _L^(51/;M"') with 
4>{l) — 0. Furthermore 



Ostp 



2\dsV 



k\ 



and therefore 



Moreover 9s0 e L^(r2i;]R^) with (9s0(O == 0. Finally, 9^80 e L^(ri;;IR^) can be concluded from 
dssVk e i2(r2,;M3)^ dsip,^ e L^i^i), dsVk £ L°°(rii;M3) and V' G i°°(f^z) using the chain rule. This 
shows e iJg J (r^j ) . D 

Remark 5. Note that e'^^^ is not injective on Hq i{flf,M.'^). Nevertheless we denote the map- 
ping ip i-^ cj) in (3.7) by (e'j,, ]^)~^, because e'^ , j^(e'j. , j^)^-'^ is the identity on iJQ;(rii). Of course, 
i^k+i)^^ ^ -^(^0 /(^Os ^0 i(^'' K'^)) &?/ ^^e inverse mapping theorem (applied in the proper quotient 
space setting). 

Theorem 6 (Existence of discrete solution). For k — 1, ...,N —1 let r^+i be the minimizer of Jk+i 
on Kk+i: provided in Theorem^ and 

^k+i :- -4+i(r/c+i)(4+i)"'- (3.8) 

Then \k+i G H^^{n.i) (see Remark\&^, and (rfe+i, Afc+i) are solving weakly the discrete fiber system 
( [33] ), i.e., 

^ (Dfc+irfe+i, (/>)2,2(n,) = -Hi^(^Q^){dsYk ■ dscf), >^k+i)H-^(n,) 

-b{dssrk+i,dss(f))L^(n,) + (ffe+i, 0)l2(o,) (3.9a) 

dsivk+i - rfe) • dsTk = 0. (3.9b) 

for all test functions G Hq i{fli;M.^). 



Proof. By definition, all elements from Kk-^-i fulfill (3.9b I. Furthermore, since r^+i minimizes Jk+i 
on Kk+i and 

efc+i[0 - (el.+i)-i[2a,r;.. • 9,0]] = for all G i/o'/^i;^^), 

we have 

= 4+,{vk+i)[4>-{e'k+i)-'[2dsrk ■ 9,0]] 

= 2(w(Dfe_j,irfe+i,0)i2(o,) +5(9s,rfe+i,as,0)z,2(f2,) - (ffc+i, 0)l2(o,) 

+ Hl,{n,){dsrk ■ dsct>, Xk+i)H-^ino) ^o"^ ^11 '^ ^ ^o,i(^i;IR')- 

D 

3.3. Stability estimates. In the following we provide stability estimates for the discrete solution. 
For function spaces 51(5*1) and 32(82) on sets Si and 5*1, respectively, we define as usual i?i(5i) (X) 
52(^2) := span{/i/2l/i € Bi(5i),/2 G -82(^2)}, where (/i/2)(si,S2) := /i(si)/2(s2), Si G Su 
S2 G ^2 (algebraic tensor product). The Sobolev space _ff^'^ (fi; K^) on fJ := 17/ x VIt then is defined as 
the completion oi H'^{ni;M.^)^H^{il.T) w.r.t. the metric associated to its inner product, see e.g. |T7l 
Chap. II. 4] (i.e., it is the Sobolev space of functions on fl which are twice weakly differentiable in 
the first variable and once weakly differentiable in the second variable and square integrable on il 
together with their derivatives). Correspondingly, we set H^l'^^{n), mi,mT G [0,C!o), to be the 
completion of H^li^i) ® if"^(rjT) and use the notation F-™'^-™r(^) := (F™/;^^(f]))'. In the 
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case 771/ = rriT we suppress the index rriT. For functions h defined on U.^ we use the following 
notation for discrete derivatives: 

(Dfe)fc ^ (Di/i)fe := ^'- ~ ^^-\ (D"/^)fc :- (D(D"-i/^)-)fe, fc = n,...,iV. 

T 

Proposition 7 (Stability estimates for r'^). Let rfe_|_i & V be as in Theorem^ k — 1, ...,N — 1, 
and let r"^ G iJ^'^(r2;M'^) be the corresponding linear interpolation. Then there exists < K < cx), 
independent of N ^ N (or, equivalently, the time discretization r > j, such that 

max ||(Dr^)fc||i2(j2,) < iC, max \\dssYk\\L'^(ni) < K, and ||r^||^2,i(j^) < JsT. (3.10) 

Proof. Since r^+i e V^ C iJ^(fi;;M'^) for all fc = 1, ...,iV — 1 and r'^(s) is piecewise linear for all 
s e rJi, we have r^ € ij2.i(r2;M3-) -^g j^j^p^^^ ^j^^t 

w ((D^r'')fe+i,(/>)i2(j2,) = - H^ni){dsYk ■ ds(}), ^k+i) h-^{q.i) 

- b{dssrk+i,dss(t>)L^{ni) + itk+i,4>)L^{ni), 

for all cj) € H^ i{ili;M.^). Note that the first summand on the right-hand side is discretized in 
an explicit way. Since (D^r'^)fc+i — ((Dr'^)i._|_i — (Dr'^)fc)/T, the special choice 4> = r^+i — r^ e 
i7^;(fiz;K3) results in 

uj (||(Dr-),+i||^.(,,,) - ||(Dr-)fe||^.(,,^) + \\r{B\-)k+i\\h^n,)) 

= -b [Wdssrk+iWh^n,) ^ \\dssrk\\l2(^ni) + \\dssirk+i ~ rfc)llL2(n,)) + 2(ffe+i, r^+i - Vk)L-^(a{) 

by applying the identity 2(a — &)a — a^ — 6^ + (a — &)^ and the functional constraint e^+i (rfe+i) = 0. 
Hence, we obtain 

w||(Dr^)fc+i||^2(o,)+fe||9^srfe+i||i2(f2,) 

<^||(Dr-)fc||i2(f,,)+6||a,.rfe||i2(a,)+2r(f,+i,(Dr-)fc+i)i2(f,,). 
Summing upfc = l,...,M— 1<A^— 1 gives the following crucial relation 

^||(Dr-)M|li2(,,,)+5||9,,rM||i2(,,,) <2r 5](ffc+i,(Dr-)fe+i)i2(o,), (3.11) 

fc=i 

(note that (Dr'^)i ~ 0, as well as, dss^i — 0). We estimate the scalar product on the right-hand 
side by Cauchy-Schwarz and Young's inequality, i.e. lab < a'^ + b^, and find 

(N-l M-1 \ 

Ellf'^-+illi^(o,)+Ell(Dr^Willi^(o,) • 
fc=l A,-l / 

The discrete Gronwall Lemma implies 

N-l ^rj.\ 

ll(Dr^)Mlli^aM ^ ^ E \\fk+i\\l^n,)e^p(-\ . (3.12) 

Together with 

N-l 
k=l 



(3.12) yields the existence oi < Ki < oo, independent of A^ G N, such that 

„T\ ||2 



(Dr-)Af lli^fo,) < Ki. (3.13) 

e existe 

(5.«r-)M||i2(o,) < i^2. (3.14) 



Combining (3.11 1 and (3.13) gives finally the existence of < K2 < 00, independent of A^ G N, such 
that 
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The inequalities (3.131, (3.14) together with the Poincare inequahty guarantee the existence of the 
desired < iiT < oo in ( |3.10[ ), independent of A^ e N. D 

Before providing the stabihty bound for the Lagrange multipher, we need to estabhsh some 
bounds on the inverse of the hnearized constraint functional. 

Lemma 8. Let (e'^+i)^^ £ L{Hli{ni);Hli{ni;R^)) be as m Remark^ k ^ l,...,N - 1. Then 
there exists < L < oo, independent of N (zN, such that 



11(4+1) ^[g]\\H^ni)<L\\g\\H^ni), ll(4+i) ^MIl^q,) < L\\g\\L2(n^), 



(3.15) 



for all g e H^i{Qi). 

Proof. Let g £ H^i{Qi). As consequence of the Poincare inequality, Lemmap^and (3.10), we find 
constants < ii, L2, £3 < oo (independent of iV e N) such that 



l(4+i)"'[5]lk^(o,)<^i 



g dsVk 

\dsVk\- 



L^(Q.i) 



< Li {\\dsgdsrk\\L^(ni) + '?>\\gdssYk\\L^(ni) 



< L2 [\\dsrk\\ca(n;) +3\\dssrk\\mn,)) \\g\\m(n,) < LsWgWmin,)- 
Using Cauchy-Schwarz inequality and Lemma [3] we obtain 



\M+ir'mhm 



gdsTk 



ds 



m^rr)k^,i9 



2 w 



< 



< 2 



gdsJTk -rfc_i) 



;i 

' 9 dsTk _ gdsVk-i 

\ 2T\d,VkV 2r|a,rfc_i|2 
2 



3^2 fg^^^v^ 3P 2 



The estimation of the discrete derivative is a little more lengthly. Integration by parts yields 

2 



ds 



du 



du 



gdsTk- 



ria^rfcP r|9^rfc_i|2 



du 



'-^^^^ dj ^.^ ] ds 



IdsTkl 



du + 2 



Yk - Vk-i g 



T IdsTk 



ds 



gdsV 



fc-i 



„ T\dsYk\'^\dsYk-l\'^ 



ds{vk - Tfc-i) • ds{rk + Vk-i) ds 



du 



<U\\g\\m(n,)+'^ 



ds\ OsijCk +rk-i)— ^TTT^ ^ I ds 



IdsTkl'^ldsTk-ll"' 



du 



Tk -Tk-l r, , , ■, gdsVk-l 

• ds{rk + Tk-i) 



ds 



< L5 I|5|Ihi(o,) 

for some < L4, L5 < 00 (independent of A^ £ N). Here, again, we used Lemmal3]and the estimates 
provided in Proposition [7] in several steps. D 



Proposition 9 (Stability bound for Lagrange multiplier). Let Xk+i £ H~^{ni) be as in (3.8), 
k = l,...,iV— 1, and let A^ be the corresponding linear interpolation. Then X^ G _ff^^'°(f2) and 
there exists < C < 00, independent of N €N, such that 



for all geH^i{ni),he i^oVl^r). 



(3.16) 
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Proof. According to (3.8 1, the Lagrange multiplier is given by 

Xk+i = - Jfe+i(rfc+i)(e'fc+i)-\ fc = 1, . . . ,iV - 1. 
Let g e Hq i{ni) and h e HQrp{^rp). Then, due to the imposed initial and boundary conditions 



tJ 

= ^/ (5, -{Xk+i-Xk) + Xk) h{t)dt 



Af-l 



2^ ffi,(a,)(5'Afc+i)i/-i(n,)'ifc+i ^ Hl,{ni){9^Xk}H-i(n, 



h 



(-1) 



k=l 



/ \ 1 /■''' + 1 



N-l 



^T.K:i^oi9,Xk)H-Hn.){^'h^-'^)k+i-HUno{9^X^)H-HnA^h(-'^)M-' (3-17) 



fc=2 



where h^ ^^ is the primitive function of h^ ^+^) with /i' ^^(T) = 0, j = 1, 2, /i^^^ = /i. Furthermore, 

Hl^{ni){9-,Xk+i)H-\i-i,) = ^•4+i(rfe+i)(4+i)" [5] 

-(ffc+i'(4+irM5])^.(^,)). (3.18) 

Using (3.17) and (3.18), now we estimate ^1 {n){g ® h, X'^)H-i(n) term by term. First we consider 



AT-l 

•E 

fc=2 



rfc - 2rfc_i - rfc_2 . , x-i 



(4)-M5]) (D^/^(-^^).+, 



< 



N--2 



r^ (DO„(D((e')-)^),,j5]) ,,^^ fD^/.( 



(-2)^ 



+ 



fc=2 
JV-2 



r5^((DO,,(e^+i)-M 



fc=2 



'L2(n, 



L2(0,) 



'fc+i 



fc+2 



+ ((Dr^)A._i,(e^_i)-Mg])^.... (d2/.(-2)) 



'TV 



JV^l 



<KL\\g\\Hiin,)rJ2 (|(D' 



2;, (-2) 



/l( 



'fe+1 



/£ = 2 



(D^'^^-^O.,! 



< 5VTA'L||5||/fi(o,)||^||ffi(Or) 
where we used several times (3.10) and (3.15). Since /i(^^^(T) = h{T) = 0, we have 



r^r — 2rAr_i — rAr_2 / , N_ir , 
-^ AeN) [9] 



L^n, 



{Bh<^- 



' N 



< CihW H^ni)\\h\\H^QT) 



for some < Ci < 00, independent of iV e N. Using (3.10), (3.15) and the continuity of f, a 



derivation of an appropriate bound for the remaining four terms is straight forward. 



D 
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3.4. Convergence to a weak solution. 

Theorem 10. There exists a sequence of discretizations (Tn)net-h r G H^'^(fl) and A G H^"{n) 
such that 



limr^"=r in C°{[0,T];L^{ni;R^)), 

limr^"=r weakly m ij2'^(17; M^), 

lim r'^" =r strongly in L^(il;M^), 

lim A^" — A strongly in H^^{n), 



(3.19a) 
(3.19b) 
(3.19c) 
(3.19d) 



for all 3/2 < /3 < oo. Furthermore, (r. A) are weakly solving (2.1 ), i.e., 

-w {dtr, dt(j))L^(n) = ~//if ^(o)(^si" ' c>s(t>, ^)ff-i.o(o) ^ ^ (^^si", dss(f>)mn} + (f , (f^)L^n) (3.20a) 



(|9.r|2,ai, 



iL^in) 



0, 



(3.20b) 



/or aZ; e i/3j(fi,;M3) ^ HIj,{VLt) and all (f) € C^i^l). Furthermore, for all < -f < 1/2 

lim r^" (t) = r(t) in C"''f{[0,l];R^), 

n— >-oo 

lim asr^"(t) = a^r(t) m C"^'^{[0,l];R^), 

n— fC30 

/or aZ^ t (^ D, where D C [0,T] is countable (in the following we are choosing D C [0,r] dense and 
let 0,T Cz D). Moreover, r has a (unique) continuous version (denoted by the same symbol) and 



r{l,t) = 0, dsr{l,t) 



r(s, 0) == (; - s)eg for all {s, t) G [0, 1] x [0, T], 



ana even 



\dsv{s,t)\^ = 1 for a.e. {s,t) € [0,^] x [0,r]. 



(3.21) 



Remark 11. (i) Theorem 10 provides existence of a weak solution to the nonlinear fourth order 

partial differential system with algebraic constraint given in (2.1 1. Taking into account the regularity 

of the weak solution, the Dirichlet and Neumann boundary conditions in ( 2.2 I are fulfilled as far as 

possible in this situation. 

(ii) Under the assumption of uniqueness of a weak solution as in Theorem 10 one has convergence 

even for all discretizations (T„)„gN- 

(Hi) The pairing j^i.o ^^JOsr ■ ds(f>, X) if-i,o(Q) in (3.20a| has to be understood in the sense of a 

uniquely closable bilinear form. I.e., if {fn)neN converges weakly to f in Hf-^'ij,{D,), {Fn)nefi *■* o 



sequence in H ' (J7) which converges strongly to F in H (fJ), lim„_>.oo ^lo ,^Jf„, Fn) fj-i,o 
exists and f — or F ^ 0, then lim„^oo ffi-." (mifn, Fn) ^^i,o 



m 



'm 



Proof. From the first two estimates in ( |3.10 ) we can conclude that r^ is uniformly (in N G 
Lipschitz continuous in C°([0, T]; ^^(f^,; M^)) and 

K(t) 1 1 e [0,r]} C {v e H^{ni;R'') \ \\dssv\\mn,) < K}, 



which is a relative compact subset of L^(ri;;M^). Thus, there exists sequence of discretizations 
(T„)„gN and r G C°i[0,T]; L^{ni;R^)) such that r^" converges to r in C°{[0,T]; L^{ni;R^)) for 
n -^ oo. 

The third estimate in (3.10 1 gives the existence of a subsequence {Tn)neN (denoted the same) and 
f e H'^'^{n;R^) such that r'^" converges weakly to f in 7J^^^(f7;M"^) for n — >■ oo. Since convergence 
in C°([0,T]; L^(ri;;K^)) implies strong convergence in L'^{n;R^) as well as weak convergence in 
i/^'^(51;M'^) implies weak convergence in L'^{n;R^), we have f ~ r. In particular, this shows 
KT9ah-K]M. 



From ( |3.16| together with the fact that the embedding H^^ (Qa) C H^i^a), a £ {I, T}, is Hilbert- 
Schmidt for all 3/2 < /3i < oo, we obtain by the kernel theorem, see e.g. [4} Chap. 1,§2.3], that 
A"^ is uniformly (in iV G N) bounded in H~^'-{il.) for all 3/2 < /3i < oo. Since the embedding 
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H^^^{Q) C H~^^~^^{Q) is compact for all < (32 < oo, there exists a subsequence (T„)„gpf and 
A e H~^[^) such that A"^" converges strongly to A in H^^{il) as n — > oo for all 3/2 < (3 < oo. 

Multiplying the linear interpolation of (3.9a) with a time-dependent test function and integrating 
w.r.t. time yields 

-~b{d,sr\dss<l>)mn) + {i\<f>)LHn), (3.22) 

for all cf) e Hq i{ili;R^) (^ H^ j,{V,t) (because Aq = Ai = 0, on [— r, 0] we can assign to the function 
dsr'^{- — t) any value, for simplicity we choose zero). Since dss^'^'^ converges weakly to dgs^ and f^" 
converges strongly to f, both in L^(ri;K'^), as n — > oo, we have 



lim (assr^",9ss0)L2(o) = (^^^r, 9ss0)l2(o), lim (F",0)L2(n) = (f, 0)l2(o), 



(3.23) 



for all (/)£ Hli{Vli 



H^rpiflx). Furthermore, integration by parts yields 



((dV)^0) 



^-1 /-tfc+i A 



L^{n) 



E 



A;=0 •"'^' 



t-tk 



(D2r-)fe+i - (D^r-),. 



(D^r-),. 



(pdsdt 



k=0 ''^ 



tk Jo 
t-tk 



T 
N-2 ,t, 



(Dr^), - (Dr^)fc_i + (Dr^)fe_i • 0dsdi 



E 



fe=0 "*'» 



t-i*. 



tjv-i "'0 





1 ft-tN_ 

T 



((Dr-)fc+i - (Dr-)fc) + (Dr^) 



0(-+t)-0 



dsdt 



ftN-l pi 



Jti JO 



-((Dr^)w-(Dr^)Ar_i)+(Dr^)Ar_i) -cfydsdi 

t)-0 



/ ((Dr^)Ar-(Dr^)Ar_i) •— / / cjidudtds + {Dr^JN^i ■ - (j)dtds 

Jo ^ / r Jt„_j J. Jo T Jt^_i 



0(. 



ds dt. 



(3.24) 



By (3.10) together with the boundary conditions imposed on cf) now from ( 3.24[ ) it follows 

Jim ((D^r-)-,0),,(,^ -^^'^^fr fo ^°'"^" • ^^^^^^^^^^^ = -(5*r,a.0),,(,, 

(3.25) 
where in the last step we used that dtv'^^ converges weakly to dtr, {4>{- + t„) — 4>)/t„ converges 
strongly to dt4> and dtv'^" — (Pr"^" ) " converges strongly to 0, all in L^(J7;M'^) as n — > oo. 
Combining ( |3.22| ), ( |3.23[ ), ( |3.25| ) we obtain 

Jl™^ffi-o^(0)(i9sr^"(- -r„) •9s0,A^")^_i.o(f2) = a; (a^r, 0)^.,^^-^^ - b{dssr,dss(p)mn) 

+ (f,0)L^(o), (3.26) 

for aU e i/2^(fi;;M3) ^ ijQi,j,(f7T). Now we restrict ourself to e H^i{ni;R^) <g) H^j^inr). 
Since dgr'^" converges weakly to dgV in Hf^'^ j,{fl;MP) and 9s0, 9ss0 are bounded functions, also 
dgr'^" (• — T„) • 9s0 converges weakly to 9^1" • 9^0 in iJg'j j.{fl) as n ^ oo. Furthermore, A'^" converges 
strongly to A, e.g. in H~^{^), as n — ^ oo. Thus 

„l™^<-,%(o)(^^'^^"(- - ^n) • ^s'/'' ^^">H-i.o(n) = Hl-l^(n){dsv ■ 9.0, A)^-i,o(s^), (3.27) 

for all g i?Q ;(57i;]R'^) (g) i/p 2^(577-) in the sense of a uniquely closable bilinear form, see Remark 
TTJiii). Hence, '([3J20a| follows from ( |3.26| ) together with p.27[ ). 
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Now, using (3.9b I, we obtain for all <j) £ C^i^) 



\i\9s 



„-r|2 



'L2(o)| 



A^-1 rt 



E 



fc=0 " '^''- 

2^1 /•*.+! fWt-tk 



E 



fc=0 •"^'» 



f (^9srfc+i - dsTkj + dsVk j dt(t)ds dt 

ids^k+i - dsVkj + ds^k ) (9srfc+i - dsvA(j)dsdt 



Jv-i .t,+i ./ 



^11 f" [ *-^(\dsrk+i\'-\dsrk\')\^\dsdt 



fc=0 ^'^ 



< 



2^E/ 



"^(-' + r^ + r^)rft||0llco(O) 



(3.28) 
fc=o •'^'= 

due to Lemma [3J Because r"^" converges strongly to r and dgV'^" , dssJ^^" converge weakly to dgV, 
dssTc, respectively, in L'^{Q;R^) as n — > oo, by an integration by parts together with (3.28) we can 
conclude 



{\dsr\^dt^ 



= - lim (d st(t) ^''" , dsr''" 

n— >-oo 



IL^{Q.) 



= 



lim ((9t0 r'^",93sr'^") 



L2(n) 



(3.29) 



for all (j) e C^in), i.e., (3.20b I is shown 



Due to the second estimate in (3.101, for each t e [0,T] there exists a subsequence (T„)„gN 
(depending on i) such that 

lim r^"(t) = r(t), lim a,r^"(i) = a,r(i) both in C°'''{[0,l];R^), (3.30) 

for all < 7 < 1/2. Let D C [0,T] be countable. Then, by dropping to subsequences and taking 



the diagonal sequence, we obtain (3.30) for all t d D. Here we choose D C [0,T] dense with 



0,T e [0,T]. From this, together with the first and second estimate in (3.10), we can conclude 
that r has a (unique) continuous version on [0,^] x [0,T] (which we denote by the same symbol). 
Moreover, 

r(/,t)=0, dsr{l,t) = -eg, r(s, 0) = (/ - s)eg for all (s,t) G [0, /] x [0,r]. (3.31) 

Finally, ( |3.29| ) together with p.31[ ) implies ( |3.21| . D 



4. Numerical study 
On every time level t^ the semi-discretized fiber system (|3.1[) corresponds to a constrained mini- 



mization problem (3.6) in a Hilbert space setting. In this section we solve the optimization problem 



numerically in a finite dimensional approximation space (finite element space) by applying a pro- 
jected gradient method with Armijo step size rule. The convergence of the method is ensured by 
the choice of a minimal projection. The temporal evolution of the fiber behavior is then realized 
by a recursive solving of the optimization problems where the solution r/j is used as initial guess 
in the projected gradient method for tk+i- The numerical results coincide well with the previous 
theoretical investigations. 

4.1. Projected gradient method on finite element spaces. Based on the classical method of 
steepest descent the projected gradient method is of first order convergence. But the theoretical 
convergence result requires the choice of a minimal projection. However, in practice also other 
projection mappings might lead to satisfying results; for details see e.g. [H]. For the fiber system in 
the Hilbert space setting the projected gradient method is given by Algorithm [l2J Note that here 
and in the following we suppress the time index ^ to facilitate the readability. 
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Algorithm 12 (Projected gradient method). 

Let P : V ^ K be a minimal projection mapping onto the set K . 

Define p :V —^ V by p(v) := v — P(v — V J(v)) via the cost functional J . 

Initialize iteration counter i = 0, choose r^^' G K, set absolute and relative tolerances tola, tolr- 

While the stationarity measure satisfies \\p{r^^>)\\ij2/Qj\ > tola + tolr\\p{y^^'^')\\H^(Qi), ^'^■ 

- calculate a^^' by the projected Armijo rule [5] such that 

J(P(rW)) < J(rW), r« = r« - a(*)VJ(r«) 

- set r(*+-^' — P{Ya ) and update i 
end 

For the numerical realization we apply the projected gradient method to finite dimensional ap- 
proximation spaces. Therefore we introduce finite element spaces Hh C H^{Q,i) and Vh C V. 
Depending on the definition of Kh we propose and investigate different projections Ph. We indicate 
all quantities associated to the finite element discretization by the space index h that should not be 
mistaken for the time index k- 

Considering the space interval [0, /], we introduce an equidistant grid {sj — {j — l)h | j = 1, ..., M} 
with space step size h — 1/{M — 1). In addition, to simplify the notation we define Sq = ■si and 
sm+1 = sm- Certainly one can also think of an irregular grid with different hj, then the partition 
is assumed to satisfy h — maxj /ij — >■ as M — )■ oo. As approximation space Hh C H^{ili) we 
choose the finite element space of piecewise cubic polynomials for function and first derivative. In 
particular 

V',(s) = <( {2-3^ + e)/^, ^ = {2s-{s,+i+s,))/h, se[s,,s,+i] 

else 

hi-l-C + e+e)/8, ^ = (2s-(s,+s,_i))//i, se[s,_i,,s,] 
cj),{s) = { h{l-^-e+e)/8, 5 = (2s-(s,+i+s,))//i, se[s,,s,+i] 
else 

j — 1, ..., Af form a node basis, i.e. 

ipj{si) = Sij, ds^j{si) = 0, (/)j(sj) = 0, dscj)j{si) = (5.y for all i,j = 1, ..., A/. 

Then any function v^ e Hh C H'^{Q,i\W^) can be represented by 

M 

Vh(s) = ^ v.V-.Cs) + v;0,(s) (4.1) 

via its coefficient tuple v = (vi, ...vm, v'j^, ...v'j^j)^ G M^*'^" with Vj,v' e M". In the finite dimensional 
fiber space Vh C Hh, the degrees of freedom reduce by 2n due to the Dirichlet boundary conditions 
given for s — I that fix the coefficients vm and v^^. Identifying the function v/j with its coefficient 



tuple V, the cost functional of (3.4) can be expressed as quadratic function in terms of the coefficients, 

J,,(v) = v^Av + b'^v + c 
A=— T + 6T", b-2T( — r-fj, c = — r^Tr, r = -2rfc + rfc_i. 

The sparse symmetric matrices T and T" have a block structure of 2 x 2, 

/ ^^ *$ \ T" - f **" **" 

- (^ ($$)T $$ j ' V (**")^ **" 

where each block itself consists of A/ x M subblocks of size n x n, e.g. 

i,j = l,...,M. 



(*$)jj = / ip^(l)jds\, (^$")y = / dssi^z dss<l>j ds \ e 
Jo Jo 
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With I e M"^" identity matrix, the subblocks are diagonal matrices, but they only contain non-zero 
entries for |i — j| < 1. The notation holds for all other blocks analogously. The time dependence 
enters Jh via the temporally earlier fiber positions f and the outer forces f . 



The linear constraint e(v) = ds(y — rt) • dsTk = ioi s € Qi (3.3) yields infinitely many 
independent conditions for the finite number of unknowns in the approximation space. To allow for 
more than the trivial solution v^ — {rk)h, we must weaken the constraint and pose the condition 
only on d space points, d < 2(M — l)n, that should be equidistantly distributed in fij. The choice 
of d is a compromise between approximation quality and numerical realization; for big d the set Kh 
shrinks to a single function that obviously fulfills the constraint on whole fi; but does not allow for 
a dynamic fiber motion {rk)h = {^o)h for all k. For small d we obtain flexibility but pay the prize 
of a lower approximation quality. Prescribing the constraint at the spatial grid points (d = M — 1) 
is a first intuitive choice. It implies the following definition for Kf^: 

Kh = {v„ e Vh I (v; - r;.) • r;. = 0, j = 1, ..., Af - 1}. (4.2) 

where r represents the coefhcient tuple to {Yk)h- 



Remark 13. Alternatively to (4.2 1, one can certainly use the finite element functions ipj, (pj to 
evaluate the linear constraint not only at the grid points Sj, but also at interior cell points. Imposing 
the linear constraint for example additionally on the cell midpoints Sj + h/2 for j = 1, ...,M — 1, 
we obtain d = 2{M — 1) and 

Kh = {v„ e Vh I (v;. ~ r;.) • r;. = 0, (v,+i/2 - r.+i/a) • r,+i/2 =0, z = 1, ..., M - 1 (4.3) 

3 1 

where Vj+1/2 :== ^(^i+i ^ ^j) ~ 4(^^+1 + ^^) "■'^'^ ''j+1/2 analogously). 

The restriction \dsYk — ds^\ < t^ on the continuous solution in K (3.6 1 (to which we refer as 
T-inequality) is implicitly contained in the set Kh by an appropriate choice of d. 

The determination of the TJ^-minimal projection mapping Ph '.Vh ^ Kh reduces to the solving 



of a linear system of equations. According to (|4.2p the coefficient tuples v' e K" of v/j e Kh satisfy 



v'- — r'- + C(r') with C(r') orthogonal complement to r' . Let {zj./}, I — 1, ..., n — 1 be basis of C'(r') 
and set v' = r' + X]i!=i '^j,i^j,i with unknowns aj E M"^^. Then, the projection 

Phiyii, 

is uniquely determined by the coefficient vector w 
solution of the linear system 

/ v|/vl/" 

_ ( **" 

q - (^ ($$" Z)T 

with Mm = and ajv/ — known due to the Dirichlet boundary conditions. The basis vectors to the 
orthogonal complements are collected in a n x A/(n— l)-matrix that is blown up to Z e jjMnx a/ (n-i) 
by M-times row-wise replication. Although, the symmetric block matrix P has to be assembled and 
decomposed only once in the projected gradient method for every time step tk, the bottle-neck in 
the computation lies in the determination of the z^ ;. It slows down the performance drastically. 
The i/^-seminorm induces a vector norm in Vh. Since in finite dimensions all norms are equivalent, 



= argmin^^ 


eKh 


\yh 


- V/1 


HHfli) 












nt vector v\j 


= 


Vl, 


■■•,VAi 


,ai,... 


,0!.M 


re 


M(2"- 


-1)*^ beir 


ig the 


*$" Z 
ZT $$" Z 


/ 
















(4.4) 




(yi 


..., 


yM ^ 


0,y'i- 


r'l,.. 


Wm 


~^'m 


-0)^ 





one might think of using instead of (4.4) the minimal projection with respect to the Euclidian vector 



norm. The orthogonal projection is comparably easy and fast to compute: 

v' • r' 

A,2(y/0=v,„ with V, =y„ < = r^ + y-- Yp^r^, i = l,...,M. (4.5) 

The Dirichlet boundary conditions are consistently fulfilled. For other definitions of Kh ( Remark [Ts]) 
corresponding projection mappings can be formulated straightforward (see [7] for ( |4.3[ )). 

The numerical realization is performed with MATLAB Version 7.0.1 on a Intel Core 2 processor, 
2GIIz. The results are scaled with Sl-units to allow for a dimensionless presentation. 
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case A: fiber at t'=0.001 for various x 
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case B: fiber at t*=0.001 for various i 




1^=0.001 

T,=0.001-2 

1^=0.0012 

Xj=0.001 
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Figure 4.2. Fiber position at t* computed with different time steps r,; = 2 ' • 
10^3^ j ^ 0,1,..., 7. Acting force in case A (left): /(s) = (s- 1)^ sin(27r(l- s)) and 
in case B (right): /(s) = - exp(-10(s - O.S)^). 
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Figure 4.3. 
(|42l with d 



Length change A/ in dependence of time step for different sets Kh.: 
M - 1 is marked by o and (|4.3[) with d = 2(M - 1) by *. The cases 



A (solid line) and B (dashed-dottcd line) refer to the different applied outer forces, 
cf. Figure lO and Table [41] 



4.2. Results and discussion. The numerical studies show the expected performance of the pro- 
jected gradient method: the strongest decay of the cost functional Jh happens within the first 
iterations, then it slows down until the tolerances of the break-up criterion are finally fulfilled. 
Thereby, the projections Ph (4.4) and P^^^ (4.5) yield similar results: the deviation of the mini- 
mized function values of J^ is of order 0(10^^"). Requiring the same number of iterations, the 
explicit projection scheme with Ph;i is much faster than the implicit one with P^. As mentioned, 
the bottle-neck of the computation is the assembly of the projection matrix P. Being interested in 
the temporal evolution of the fiber and hence in an efficient computation of thousands of time levels 
via the projected gradient method we apply Ph.2 in the following simulations. 



We consider a set-up with eg = — ei and f = /e2 for different exemplary force functions /. 
The fiber parameters are Z = 1, w = lO^'^ and h — 10^^. The space discretization is chosen to 
be M = 300. The projected gradient method is initialized with the solution of the previous time 
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case A: / 


(s) = (s-l)3sin(27r(l-s)) 








efficiency 


accuracy 




# iterations (0, E) 


CPU-time [s] 


|(r--rO(f*)| 


A/ 


To 


1371.0 1371 


86 


2.0- 10 1 


1.6-10-2 


Tl 


342.0 684 


40 


1.0-10-1 


1.0 


10-2 


T2 


84.3 337 


19 


4.7-10-2 


7.3 


10-3 


Ta 


21.0 168 


11 


2.2- 10-2 


6.0 


10-3 


T4 


3.3 52 


7 


1.1-10-2 


5.3 


10-3 


T5 


3.6 116 


14 


5.3-10-3 


4.0 


10-3 


Ta 


3.8 244 


34 


2.1 10-=* 


3.3 


10-3 


T? 


4.2 542 


80 





2.7 


10-3 





case B: /(s) = - exp(-10(s - 0.5)2) 










efficiency 


accuracy 




# iterations 


(0, S) CPU-time [s] 


|(r-.-r-0(i*)| 


a; 


To 


1374.0 


1374 88 


5.3- 10-1 


3.1-10-2 


n 


356.0 


712 41 


2.6-10-1 


1.7 


10-2 


T2 


94.0 


376 21 


1.3-10-1 


1.2 


10-2 


T3 


24.0 


192 12 


5.5 10-2 


1.0 


10-2 


T4 


3.8 


60 7 


2.8-10-2 


8.6 


10-3 


T5 


4.1 


132 16 


1.2-10-2 


7.4 


10-3 


T6 


4.5 


289 36 


5.6 lO-'' 


6.3 


10-3 


T? 


5.1 


648 92 





5.2 


10-3 



A/ ^^ 



1.6 
5.5 
2.2 
9.8 
4.6 
2.2 
1.0 
5.3 



10^ 

10-3 

10-3 

lO-'^ 

lO-'^ 

lO-'^ 

lO-'^ 

10-5 



a; (4.3|) 


3.1 - 10-2 


9.9 


10-3 


3.8 


10-3 


1.7 


10-3 


8.1 


10-4 


3.9 


10-4 


1.9 


10-4 


9.7 


10-5 



Table 4.1. Influence of time step on performance of algorithm for Kh of (4.2 1, 
i.e. averaged number of iterations in projected gradient scheme 0, total number of 
iterations S, CPU-run time in seconds, approximation error in final fiber position, 



fiber elongation AL In addition AZ for K^ of (4.3) on the right 



level, the tolerances are set to be tola = 10 3^ tolr = 10 2. Accuracy and efficiency of the whole 
algorithm are strongly affected by the choice of the time step r. To get a qualitative impression, we 



define a sequence of time steps (Ti)^, ti — 2 
for t* = 10-3, see Figure O and Table [JT 



10 



-3 



0, 1, 2, ... and determine the fiber behavior 



In accordance to the theory, we clearly observe the 
convergence rate of the implicit Euler time discretization as first order; the approximation error is 
linear in r which is emphasized by the marked values for tq, t^, tq in Table [4T] Considering the 
computational effort, the averaged number of iterations required in the projected gradient method 
reduces for smaller r due to the better initial guess, whereas the total number of projected gradient 
runs increases. Thus, moderate choices of t (here: r ~ 1Q~^) yield the fewest number of iterations 
in total and hence the fastest computation. The fiber elongation A/(i*) = J \dsY'^{t*)\ds — / > is 
originated in Lemma [3] and reduces for smaller r, certainly AZ — > for r ^ 0. However, its actual 

Replacing the set of 

2(Af 



magnitude is a consequence of the definition of Kh- Replacing the set of ( |4.2[ ) by the one of (4.31 
where the linear constraint is imposed at the double number of space points {d 



1)) we 

obtain a better length conservation (smaller A/); in particular. A/ decreases here linearly in r, see 
Figure |4.3| It is worth to mention that the computational effort stays comparably the same when 



a respective projection Ph^2 for (4.3) is used. The reason lies in the same underlying finite element 
space Vh- 

We are interested in using the numerical scheme for the simulation of a fiber dynamics. Therefore, 
it is important that the elongation is restricted over time. Theoretically, 



Al{t)<tTl<TTl, te[0,T], 



(4.6) 



holds according to Lemma [3] The bound results from the r-inequality in the definition of K. In 
our finite element formulation, the inequality constraint is incorporated in K^ by the choice of d 
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fiber motion, d=M+1 
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Figure 4.4. Fiber dynamics over time for different sets K^ where the hnear 
constraint is satisfied at spatial points Gj^ j — l,..,d. From left to right: aj = 
ij-l)h and d = M-l (|X2|; aj = {j-l)h/2 and d = 2(A/-f) ^^: aj = {j-l)h/2, 
andd = 3(M-l). 



letigth change over tifne in dependence of K^^ 
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Figure 4.5. Elongation over time in dependence of the set Kh, cf. Figure [44} 



The dashed line visualizes the theoretical bound A/(i) < trl (4.6) 



(Remark 13). This implicit treatment avoids the costly and complex optimization with inequality 
constraints, but only allows for an a posteriori check of (4.6). We demonstrate the impact of K^ 



on the fiber behavior for an easy, meaningful scenario: in the stated set-up we apply an additional 
vertical force, i.e. f = /lei + /2e2 with /i and /2 of same order. Figure [T4| shows the fiber motion 
over time computed with t ^ t^ f=i 10""* for d e {M - 1, 2(M - 1), 3(M - 1)} (here: /i = -Wuj, 



1 



(4.2) yields an obviously wrong, 
= 2,3 look very 



h{s) = -10-2sin(27r(l - s))). Whereas Ku with d = M 

disturbing and even growing elongation, the results to Kh with d = i{M — 1), i 

similar and reasonable on the first glance. However, only the set Kh where the linear constraint 

is prescribed at the 3(M — 1) space points aj = {j — l)/i/3 fulfills the theoretical bound (4.6) and 

provides the desired numerical solution, Figure [43] This stays true for the longtime behavior. Note 

that the computational effort of the algorithm with "Euclidean projections" is determined by time 

and space discretization (and not by d). 

Remark 14. The commercial software tool FIDYST (Fiber Dynamics Simulation Tool r\ developed 
by the Fraunhofer Institute of Industrial Mathematics (ITWM) in Kaiserslautern, Germany for the 



simulation of technical textile production processes computes the fiber dynamics ( 2.1 ) by help of a full 



^FIDYST, Fraunhofer ITWM, Kaiserslautern, http : //www. itwm. fraunhofer .de 
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discretization: implicit Euler method in time and finite difference method of higher order in space. 
The resulting nonlinear system is then solved by a Newton method. In comparison to that approach 
our proposed scheme is computationally competitive. Moreover, we provide a strict theoretical basis 
/justification, including convergence result and error bound. 

5. Conclusion 

In the technical textile industry the dynamics of an elastic inextensible slender fiber is modeled 
by an nonlinear fourth order partial differential algebraic system. In this paper we have proposed 
a numerical scheme focusing on an accurate and efficient treatment of the algebraic constraint. 
Ongoing work deals with the extension of analysis and numerics to the stochastic partial differential 
algebraic system 15J arising for fibers immersed in turbulent air flows. Here, a stochastic force 
(source term) of a white noise type is added in the fiber equation. The challenge lies again in the 
handling of the constraint. So far, the corresponding extensible equation with additive Gaussian 
noise has been studied in [3]. 
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